rm(list = ls())
# Function to load packages
loadPkg=function(toLoad){
for(lib in toLoad){
if(! lib %in% installed.packages()[,1])
{install.packages(lib, repos='http://cran.rstudio.com/')}
suppressMessages( library(lib, character.only=TRUE))}}
# Load libraries
packs=c("maps", "rworldmap", "maptools", "sp", "spdep", "gstat", "cshapes","readstata13", "data.table",
"spatstat", "DataCombine", "countrycode",  "cshapes", "arm", "rgeos", "raster", "gpclib", "tidyr", "pwt9",
"texreg", "reshape2","ggplot2", 'foreign', 'car', 'lme4', 'dplyr', "mapdata", "rgdal","gridExtra",
"pscl", "stargazer","stringi","multiwayvcov", "lmtest","scales","stringi","haven",
"coefplot", "ggthemes", "gridExtra", "grid", "cowplot", "GGally")
loadPkg(packs)
setwd("YOURPATH/Replication")
getwd()
### I make a set of R functions for plotting. To replicate the all figures, you will
#need to install/load the R function in the source files ("postregplots_0.1.0.tar.gz")
##to make the plots with a simulation approach
library(postregplots)
#set ggplot theme
theme_set(theme_gray())
load("./Data/paperdata.RData")
data <- paperdata %>%
dplyr::filter(weaker_chal==1)
###########################################################################
############### This is the main results in the analysis (Models 1 - 6)
###########################################################################
#Model 1: base
b1 <- bayesglm(initiation ~ tgt_landneighbors2 +  backdown_terr1 +  backdown_nonterr1 +
joint_demo +  defensetgt + offensechal +  land_neighbor +
vicratio +  tmid +  tmid_sq + tmid_cub,
data = data, family = binomial(link = "logit"),
prior.scale=2.5, prior.df=1)
# model 2:
b2 <- bayesglm(initiation ~ tgt_landneighbors2 +  backdown_terr1 +  backdown_nonterr1 +
joint_demo +  pot_milicapchal +  pot_milicaptgt +  land_neighbor +
tmid +  tmid_sq + tmid_cub,
data = data, family = binomial(link = "logit"),
prior.scale=2.5, prior.df=1)
#model 3
b3 <- bayesglm(initiation ~ tgt_landneighbors2 +  backdown_terr1 +  backdown_nonterr1 +
joint_demo +  intra_wartgt +  inter_war_totaltgt + defensetgt +
offensechal +  land_neighbor +
vicratio +  tmid +  tmid_sq + tmid_cub,
data = data, family = binomial(link = "logit"),
prior.scale=2.5, prior.df=1)
# Model 4:
b4 <- bayesglm(initiation ~ tgt_landneighbors2 +  backdown_terr1 +  backdown_nonterr1 +
joint_demo +  solschangetgtdum + solschangechaldum +
offensechal +  land_neighbor + defensetgt  +
vicratio +  tmid +  tmid_sq + tmid_cub,
data = data, family = binomial(link = "logit"),
prior.scale=2.5, prior.df=1)
# Model 5
b5 <- bayesglm(terrmidl ~ tgt_landneighbors2 +  backdown_terr1 +  backdown_nonterr1 +
joint_demo +  offensechal +  land_neighbor + defensetgt  +
vicratio +  tterr + tterr2 + tterr3,
data = data, family = binomial(link = "logit"),
prior.scale=2.5, prior.df=1)
#Model 6:
b6 <- bayesglm(initiation ~ tgt_landneighbors2 +  backdown_terr1 +  backdown_nonterr1 +
joint_demo +  offensechal +  land_neighbor + defensetgt  +
weaker_chal +  major_major +  major_minor +  minor_major +
vicratio +   tmid +  tmid_sq + tmid_cub,
data = paperdata, family = binomial(link = "logit"),
prior.scale=2.5, prior.df=1)
################### Figure 2(a)###################
coef_f1 <- coef_clusterplot(ModelResults = b1,
varnames = c("tgt_landneighbors2","backdown_terr1",
"backdown_nonterr1", "joint_demo",
"defensetgt", "offensechal",
"land_neighbor","vicratio"),
data = data,
clusterid = "ddyad") +
scale_x_discrete(labels=c("Past yielding (nonterritory)", "Past yielding",
"Defensive alliance (target)", "Joint democracy",
"Land neighbor", "Offensive alliance (challenger)",
"Num. of land neighbors",
"Exp. victory prob. ratio")) +
theme(axis.text = element_text(size=14),
plot.title = element_text(hjust = 0.5,size = 30, face = "bold"),
axis.title=element_text(size=14)) +
ggtitle("")
coef_f1
ggsave("./figures/cof1.pdf", plot = coef_f1,  width = 10, height = 8)
################### Figure 2(b)###################
coef_f2 <- coef_clusterplot(ModelResults = b2,
varnames = c("tgt_landneighbors2","backdown_terr1",
"backdown_nonterr1", "joint_demo",
"pot_milicapchal", "pot_milicaptgt",
"land_neighbor"),
data = data,
clusterid = "ddyad") +
scale_x_discrete(labels=c("Past yielding (nonterritory)", "Past yielding",
"Joint democracy",
"Land neighbor", "Alliance strength (challenger)",
"Alliance strength (target)",
"Num. of land neighbors")) +
theme(axis.text = element_text(size=14),
plot.title = element_text(hjust = 0.5,size = 30, face = "bold"),
axis.title=element_text(size=14)) +
ggtitle("")
ggsave("./figures/cof2.pdf", plot = coef_f2,  width = 10, height = 8)
################### Figure 2(c)###################
coef_f3 <- coef_clusterplot(ModelResults = b3,
varnames = c("tgt_landneighbors2","backdown_terr1",
"backdown_nonterr1", "joint_demo",
"defensetgt", "offensechal",
"intra_wartgt", "inter_war_totaltgt",
"land_neighbor","vicratio"),
data = data,
clusterid = "ddyad") +
scale_x_discrete(labels=c("Past yielding (nonterritory)", "Past yielding",
"Defensive alliance (target)","Interstate war (target)",
"Civil war (target)", "Joint democracy",
"Land neighbor", "Offensive alliance (challenger)",
"Num. of land neighbors",
"Exp. victory prob. ratio")) +
theme(axis.text = element_text(size=14),
plot.title = element_text(hjust = 0.5,size = 30, face = "bold"),
axis.title=element_text(size=14)) +
ggtitle("")
ggsave("./figures/cof3.pdf", plot = coef_f3,  width = 10, height = 8)
################### Figure 2(d) ###################
coef_f4 <- coef_clusterplot(ModelResults = b4,
varnames = c("tgt_landneighbors2","backdown_terr1",
"backdown_nonterr1", "joint_demo",
"defensetgt", "offensechal",
"solschangetgtdum", "solschangechaldum",
"land_neighbor","vicratio"),
data = data,
clusterid = "ddyad") +
scale_x_discrete(labels=c("Past yielding (nonterritory)", "Past yielding",
"Defensive alliance (target)",
"Joint democracy",
"Land neighbor", "Offensive alliance (challenger)",
"SOLS change (challenger)",
"SOLS change (target)",
"Num. of land neighbors",
"Exp. victory prob. ratio")) +
theme(axis.text = element_text(size=14),
plot.title = element_text(hjust = 0.5,size = 30, face = "bold"),
axis.title=element_text(size=14)) +
ggtitle("")
ggsave("./figures/cof4.pdf", plot = coef_f4,  width = 10, height = 8)
################### Figure 2(e)###################
coef_f5 <- coef_clusterplot(ModelResults = b5,
varnames = c("tgt_landneighbors2","backdown_terr1",
"backdown_nonterr1", "joint_demo",
"defensetgt", "offensechal",
"land_neighbor","vicratio"),
data = data,
clusterid = "ddyad") +
scale_x_discrete(labels=c("Past yielding (nonterritory)", "Past yielding",
"Defensive alliance (target)", "Joint democracy",
"Land neighbor", "Offensive alliance (challenger)",
"Num. of land neighbors",
"Exp. victory prob. ratio")) +
theme(axis.text = element_text(size=14),
plot.title = element_text(hjust = 0.5,size = 30, face = "bold"),
axis.title=element_text(size=14)) +
ggtitle("")
ggsave("./figures/cof5.pdf", plot = coef_f5,  width = 10, height = 8)
################### Figure 2(f)###################
coef_f6 <- coef_clusterplot(ModelResults = b6,
varnames = c("tgt_landneighbors2","backdown_terr1",
"backdown_nonterr1", "joint_demo",
"defensetgt", "offensechal",
"land_neighbor","vicratio",
"weaker_chal", "major_major", "major_minor",
"minor_major"),
data = paperdata,
clusterid = "ddyad") +
scale_x_discrete(labels=c("Past yielding (nonterritory)", "Past yielding",
"Defensive alliance (target)", "Joint democracy",
"Land neighbor","Major--Major powers","Major--Minor powers",
"Minor--Major powers",
"Offensive alliance (challenger)",
"Num. of land neighbors",
"Exp. victory prob. ratio",
"Weaker challenger dyad")) +
theme(axis.text = element_text(size=14),
plot.title = element_text(hjust = 0.5,size = 30, face = "bold"),
axis.title=element_text(size=14)) +
ggtitle("")
ggsave("./figures/cof6.pdf", plot = coef_f6,  width = 10, height = 8)
b0 <- bayesglm(initiation ~ joint_demo +  defensetgt + offensechal +
vicratio +  tmid +  tmid_sq + tmid_cub,
data = data, family = binomial(link = "logit"),
prior.scale=2.5, prior.df=1)
roc <- plot_roc(ModelResults = list(b0, b1),
linetypes = c("solid", "dotted"), interval=0.2)
ggsave("./TEX/figures/roc.png", plot = roc)
ggsave("./figures/roc.png", plot = roc)
###################### Figure 3a-3b: Predicted effects of past concession
## and future stake on territorial dispute initiation by weaker states
#based on Model 1
sim_b1 <- sim_cluster(ModelResults = b1, n.sim = 1000, data = data,
clusterid = "ddyad")
###Figure 3(a)
##set simulation
ModelResults = b1
val1 =0; val2 = 20; intervals = 1
varname2_val = seq(val1, val2, by = intervals)
varname2 = "tgt_landneighbors2"
df <- array(NA, c(length(varname2_val), 4))
for (i in 1:length(varname2_val)){
X1 <- model.matrix(ModelResults)
X1[, varname2] = varname2_val[i]
fd = apply(apply(X1, 1, function (x) plogis(sim_b1 %*% x)), 1, mean)
df[i, 1] <- varname2_val[i]
df[i, 2] <- mean(fd)
df[i, 3:4] <- quantile(fd, probs = c(.05,.95))
}
colnames(df) <- c("X", "mean", "lo", "hi")
df <- as.data.frame(df)
p_b1 = ggplot(df, aes(x=X, y=mean)) +
theme_gray() +
geom_ribbon(aes(ymin = lo, ymax = hi), alpha=.2, color = NA) +
geom_line(aes(y = mean, x = X), size = 0.5, alpha = 0.75) +
scale_x_continuous(name = "Num. of land neighbors",
breaks = seq(val1, val2, by = 1)) +
theme(axis.text = element_text(size=12),
plot.title = element_text(hjust = 0.5,size = 30, face = "bold"),
axis.title=element_text(size=12)) +
ylab("Pr(initiation)")
p_b1
rm(list = ls())
